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Abstract 

This article discusses computational techniques for simulating natural con- 
vection in three-dimensional domains using finite element methods with tetra- 
hedral elements. These techniques form a new numerical procedure for this 
kind of problems. In this procedure, the treatment of advection by a wave 
equation approach is extended to three-dimensional unstructured meshes with 
tetrahedra. 

Numerical results of natural convection of an incompressible Newtonian fluid 
in a cubical enclosure at Rayleigh numbers in the range 10^ to 10^ are obtained 
and they are in good agreement with those in literature obtained by other 
methods. 

key words Three-dimensional domains, Finite element methods, Tetrahedral 
elements. Natural convection, Incompressible fluids. 

§1 Introduction 

Simulation of natural convection flows in three-dimensional geometries has been an 
area of active research in recent years. In the past decade, most researchers who 
performed calculations in three-dimensional geometries were hindered from applying 
sufficient resolutions, by limitations on computer storage. For example, Mallinson 
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and De Vahl Davis in 1977] used a very coarse mesh with up to 15^ nodes, Pepper 
D.W. in 12 1987] apphed only 33 x 17x 9 nodes. More recently, with the availability of 
more powerful computers, researchers are now able to perform calculations on meshes 
with better resolutions. Le Peutrec and Lauriat 1990] used mesh with up to 41^ 
nodes, Fusegi et al. in [H 1991] used meshes with up to 62'^ nodes. Janssen et al.'m 
13 1993] reported results with 120^ nodes mesh, but they generated the results by 
symmetry, after performing actual simulation with only one-fourths of this number. 
Natural convection is governed by a coupled system of Navier-Stokes equations and 
energy equations. 

The objective of this paper is to present a finite element method for simulating natural 
convection of an incompressible fluid in three-dimensional geometries using tetrahe- 
dral elements with unstructured mesh. An operator-splitting scheme of Marchuk- 
Yanenko is applied to split the coupled system into three sub-problems namely, the 
pressure, transport and diffusion sub-problems. This decouples the difficulties usu- 
ally associated with non-linearity and incompressibility constraint. The pressure and 
diffusion sub-problems are time discretized by backward-Euler-type method. The 
non-linear advection is treated by a wave equation approach. Space discretization is 
achieved by a finite element method where pressure, velocity and temperature are 
approximated by continuous piecewise-linear polynomials on meshes consisting of 4- 
node tetrahedral elements. The mesh for velocity and temperature is twice finer than 
the pressure mesh so that the inf-sup condition is satisfied. A systematic method of 
constructing these velocity-pressure meshes, such that each pressure tetrahedral ele- 
ment is a macro-element consisting of eight sub-tetrahedra for velocity, is discussed 
in this article. We extend the two-dimensional method for constructing a pressure 
macro-element, by connecting edge mid-points, to three dimensional meshes. The 
numerical procedure presented in this article also extends the treatment of advection 
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by a wave equation approach in [H], [7] to three-dimensions while combining other 
different numerical techniques and thus forming a new, efficient, solution procedure 
suitable for simulating motion of an incompressible fluid in three-dimensional geome- 
tries with unstructured meshes. 

Results obtained for the numerical example of natural convection of air in a cubi- 
cal box, illustrate the accuracy and reliability of this new procedure. The three- 
dimensional results also validate the usual assumptions in two-dimensional simula- 
tions and elucidates three-dimensional effects on this flow phenomenon. 

§2 Governing Equations for Natural Convection 

We consider natural convection of an incompressible viscous Newtonian fluid enclosed 
in a three-dimensional rectangular domain, Q C R^, with boundary denoted by F. 
The geometry and coordinate system for the enclosure are shown in Figure ^ The 
natural convection is induced by the non-zero temperature gradient between the two 
vertical surfaces, F/ (at a; = 0) and F^ (at x = L^). The remaining four surfaces, 
(F \ F/ U Tr) are assumed to be perfectly thermally insulated. 




Figure 1: Geometry and coordinate system for cubical enclosures. 



(0<a;<L^; < y < Ly ; 0<z< L,). 
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With the Boussinesq approximation, the vector form of the dimensionless governing 
equations in a finite time interval (0, t'^) are: 



du 

'dt 

de 
di 



+ (u ■ V)u - Pr Au + Vp = RaPrOj inQx (0, r ), 



+ (u-V)^ - = 



in Q X (o,r), 



(2) 



(1) 



Vu 







in nx (o,r), 



(3) 



where , 



Ra 



g(3(n-T,)Ll 



, the Rayleigh number, 



(4) 



Pr 



— , the Prandtl number, 

a 



(5) 



\i{x, y, z, t) — {ui}^^^ — {ux, Uy, Uz) is the flow velocity, 

t is elapsed time, (u ■ V)u is a symbolic notation for the non-linear vector term 



9{x, y, z, t) is the temperature, 
p{x, y, z, t) is the pressure , 

P is the coefficient of thermal expansion of the fluid, 
a is the coefficient of thermal diffiisivity of the fluid, 
1/ is the kinematic viscosity coefficient of the fluid, 
g is the gravitational acceleration. 

This set of dimensionless equations is subject to the following initial and boundary 
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conditions: 



(6) 



lu(x, y, z, 0) = in fl, 

6{x, y, z, 0) = in Q, 
u{x, y, z,t) = on r X (0,t^), 
e{x, y, z,t) = 1 on r, X (0,t^), 
e{x, y, z,t) =0 on T,. x (0,t^), 
^ ||(a;, y,z,t) = on T \ {Ti U T,) x (0, t^). 

In order to obtain the dimensionless equations, the distance between the colder and 
hotter surfaces, L^., has been chosen as the reference length and the scale factors for 
velocity, time and pressure are chosen as, a/L^, ja, p^a^ jL^ respectively. The 
dimensionless temperature is defined as 6* = ^ . 

Thus, the three-dimensional rectangular model domain has dimensions Ay^ A^ 
where A^^. is the aspect ratio in the Xi direction. 

§3 Time Discretization by Marchuk-Yanenko-Type Operator Splitting 
Method of Problem Equations(|H) - © 

Let At > be the time step and t"" = nAt. At every time interval [t", t"'*'"^], the 
Marchuk-Yanenko-type operator-splitting method involves a sequence of computa- 
tions as follows: 

(I) = uo ; ^° = ^0, (7) 
then, for n>0, u",^" given, we compute (e^+^Z^^ u"+i/^p"+i) , (^"+2/3^ ^"+2/3^ 

(Qn+l n+l\ follows: 
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(II) Solve the pressure sub-problems: 



Qn+l/3 ^ Qn -^^ 
n+1/3 



U 



U 



Vp"+^ in 



^n+i/s = on r, 
V ■ u"+i/3 = in 1], 

(III) then solve the transport sub-problems 

^ II + u"+i/3 . ^ in ^ X (r, 
^ + u"+i/3 . ^ in X (r, t 

^ u(r) = u"+i/^ ^(r) = 



+1^ 



71 -tn+l^ 



(9) 



Set u"+2/3 = _ ^(^n+l)_ ^^q) 

(IV) Finally, solve the diffusion sub-problems: 

' nn+l /in+2/3 

2 ^1 = inl], 



in+l 



1 on Y f 



r+i = on r^, 

1^"^' = on r\(r,ur,), 

„n+l ,,'"+2/3 

^ ^ Pr Au"+i = Ra Pr ^"+1 j in fi, 

u"+i = on r. 



(11^ 



§4 On the Finite Element Approximation of sub-problems (jSI) - (jlip 

The entire domain Q = Q^, constituting the computational domain, is discretized 
into a finite set, 7^, of tetrahedra inside which velocity, pressure and temperature are 
continuous and the collection of the tetrahedra satisfies the following properties: 

1. Th = Th G Qh for all Th G %. 



2. %^ {Th} is finite. 

3. For Ti , T2 e 7/i, if interior (Ti) ^ interior(T2), then only one of the following 
is possible: 

• Ti n TIj = a vertex common to Ti and T2 or, 

• Ti n T2 = a common edge or face of Ti and T2 or, 

• Ti n ^2 = 0- 

4. U =^^/.- 

The finite element mesh for pressure is twice coarser than the mesh for velocity. The 
mesh for temperature is the same as that for velocity. Let 7^ be the finite collection of 
the tetrahedra for velocity and 72/i, a similar collection for pressure. Piecewise-linear 
approximation is employed for all variables including pressure. Hence, the vertices 
of the tetrahedra in Th and form the nodes for the finite element meshes for 
velocity,temperature and pressure respectively. 

A new method of discretizing Vth and Vt2h into 4-node tetrahedral elements, such 
that the pressure tetrahedra are macro-elements consisting of eight sub-tetrahedra, 
is presented in this report. This is achieved by first discretizing the domain into a 
finite set of 8-cornered brick-like macros, each of which is then split into two prisms 
along a vertical mid-plane through a diagonal line on the top face and through the 
center of gravity of the brick-like macro. Each of the two prims is then divided into 3 
tetrahedra in a very unique way. is constructed first and, as in two-dimensional 
cases, Th is constructed from by connecting the edge midpoints on the faces and 
on any resulting vertical mid-plane. A detailed explanation of the steps involved 
in this tetrahedralization is given in the appendix. We however make the following 
remarks here: 
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§5 Remarks 

1. The six generic tetrahedral elements have different shapes but the same volume 
which is one-sixth of the volume of the particular brick-like macro containing 
them. 

2. In order for each pressure tetrahedral element to properly contain exactly eight 
sub-tetrahedra, after the edge-midpoints have been connected, it is necessary 
that, one set of three tetrahedra on one prism must be a reflection of the second 
set of three tetrahedra on the other prism about the vertical plane through the 
diagonal of the brick-like macro. 

3. Although, a method, for splitting a brick-like macro into six tetrahedra was 
given, by Zienkiewcz, in jSj, the discussion did not include the situation of 
multi-level grids as in this present case. 

§6 Discrete sub-problems and Weak Formulations 

Weak formulation of each set of sub-problems, determined by the operator-splitting, 
are obtained using the following fundamental discrete spaces: 

Vh = {vH IvheC" (n,) , v,,\t e Pi, V T G %} , (12) 

Voh = {vh I Vh e Vh, Vhlr = 0} , (13) 

Tloh = {(j)h\(j) e Vh, Mr.urr = 0} , (14) 

Ph = \qh I g/. e C° {n2h) , qhW eVi,W e T^h, f qtdn = 1 . (15) 

In eg nations (fT^ - (fT^ . Vi is the space of polynomials in three variables of degree 
< 1. The discrete approximation associated to the flnite element spaces described 
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above, for the weak formulation of the pressure sub-problems is: 
Find K+'/^ p^+i} G X such that, 

1 f .."+1/3.. ro _ /" ^n+l9Vh,r^ , 1 



At 



f uT^'^VhdQ= I pl+^^dQ + ^ l uyan ^VheVoh, (16) 

for each component of \ih = {ui}^^^ , 

g;,V-u;;+'/'dl] = 0, yqhePh, (17) 



u;;+i/3 = o, on r. (18) 

The transport sub-problems combined with the wave equation approach described in 
jU] and j7] give a set of semi-discrete sub-problems with weak formulation given as: 
Find Oh e Vh and Uh G (K)/^)^ Vt G [r, such that, 

/ f ^) <l>hdn + f (u^^/^ ■ V0,) (u^^/^ • v^,) dn = o, V 0, G y,, (i9) 

/ (^) ■^hdQ+ [ (u^^/^ ■ Vw,) (u^^/=^ ■ Vu,) rfl] = 0, V w, G (\4)' , 

(20) 

^/.(i") = C^', Mn = <^"\ (21) 

{1 on F/., 
(22) 
on Fr 

The solution of a wave-like equation such as eq. (fT^ and eq. fpU]) has been described in 
in] and [7j for uniformly structured meshes. It involves time discretizing the equation 
by a second-order finite difference scheme with an initialization step consisting of 
solution of full discrete version of eq.©- It is noteworthy to mention that since the 
discrete wave-like equations are explicit there is no need to store any square matrix 
thereby, conserving computer memory. However, a local time step At/Q has to be 
chosen with integer Q sufficiently large so that the CFL condition is not violated. 
The set eq. (fTT|) . of diffusion subproblems is approximated by the following discrete 
sub-problems: 
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Find ei^^ e Vh and G such that, 

kdn + [ \/ei+'-V(f)hdn = o, \/ (f)^ e no,., (23) 



nn+l _ nn+2/3 
h h \ A. \ I V7/3"+l 



At 

, on Fr, 

= (24) 



1 on F 



h 



At 



■Whdn +Pr [ Vul^^ -Vwhdn 

= RaPr[ ei'^rwhdn, e {Vh)\ (25) 

§7 Solution Strategy for the Pressure Subproblems (jl6|) - (jlSp 

Let X = (x, 2;) and {ciJm(x)}^^^ be the vector basis for Vh such that, 

{1 if = Xm 
(26) 
if x^ 7^ Xm 

for each node on the velocity mesh. After applying Galerkin method on the discrete 
weak formulation we get the following discrete sub-problem at each node x^ on the 
velocity mesh: 

Find K+'/^K+l} G (K)0' X n such that, 

[ ^C'uJmdn= I pl+^^dn+^f w>^dn, (l<z<3) (27) 

u;;+^/' = 0, on F, (28) 
UkV ■ Uh^^^^ dQ = 0, at each node x^ on the pressure mesh. (29) 

This sub-problems is solved by a preconditioned conjugate gradient (PCG) algorithm 
described in jH] and Proper evaluation of some integrals in the PCG algorithm 
is very crucial to the overall performance of this numerical technique. These include 
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integrals of the forms f p^^^ dQ and / cu/t V ■ u^^ dQ which involve product 
of discrete functions over coarse and fine meshes. The following is a summary of the 
techniques applied to these integrals: 



1. All integrals of the form / fgdVt are approximated by the trapezoidal method 

Jn 

globally on the element domain. 



2. All integrals of the form / p^—^dVL are computed on the fine velocity mesh 
element-by-element. On each pressure element, the function ph is interpolated 
linearly along element edges. Also, over each velocity element, ph is approxi- 
mated by the average of its nodal values on the vertices. 



3. The integrals of the form / cj^ V ■ u/^ dVt are computed over the coarse pressure 
mesh element-by-element. In each pressure tetrahedral element, the nodal val- 
ues of the coarse-mesh basis function uj^ on the vertices of each of the 8 included 
velocity elements are obtained by linear interpolation along the edges. V ■ u/j 
is piecewise-constant over each velocity element since \ih is approximated by a 
piecewise-linear function. 

The PCG steps include the solution of a Neumann problem with solution belonging 
to the space of functions with mean-value zero. Since this problem is solved on a 
relatively coarser mesh, with banded storage, it is expected that memory requirement 
will be a manageable size on most modern computers. Thus the Neumann problem is 
solved by direct method after cholesky factorization as suggested by Glowinski in j^l 
page 267]. However, since its solution has mean- value zero the following steps must 
be performed together with the direct method: 

1. Set one of the unknowns to zero and delete the corresponding row and column; 
The NxN linear system in (say) will reduce to (N — 1) x (N — 1) linear system 
in 0', 
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Solve the reduced linear system by direct method, for 0', 



3 



compute mean value of (j)' = — — / (p' dVt, 

measiil) I 
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§8 Solution Strategy for the Diffusion Subproblems (j23j) and (j25j) 

The integrals in eq. ()23|) and eq. ()25p are assembled over all the fine-mesh tetrahedral 
elements. With proper ordering of nodes, the diffusion sub-problems for temperature 
and velocity components result in linear systems of nodal values of the form 



where A is a symmetric positive definite, banded, sparse matrix. The bandwidth 
of A, however grows very rapidly as the resolution is increased, so that memory 
requirement becomes prohibitively large even on supercomputers, despite banded 
storage. Thus it is more practicable to solve these linear systems for temperature 
and velocity components by iterative methods. A careful observation reveals that 
with piecewise-linear approximations, A has only 15 non-zero diagonals in each case. 
Since it is also symmetric, a substantial amount of memory is freed by storing only 
the main diagonal and the 7 non-zero upper (or lower) diagonals. 
The associated Dirichlet boundary condition should be enforced in a manner that pre- 
serves the symmetry and sparse nature of A. A method of achieving this is discussed 
by Stasa in [TU[ pp 59-61]. The diffusion linear systems of the form eq. ljHnj) are solved 
by conjugate gradient algorithm of Hestenes and Steifels (CGHS) given for example 
in Plj. This algorithm involves matrix- vector multiplications which are performed 
within the bandwidth and in a manner that prevents ^^fill-ins" . Only 7 multiplication 



Ax = f 



(30) 
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operations are required per row. To accelerate and ensure convergence, precondition- 
ing is usually necessary. A way to achieve this is to ensure that Gerschgorin disks 
for the iterates are concentric by performing symmetric scaling, where the main di- 
agonal elements are scaled to unity before commencing the iteration process. The 
linear system for temperature is solved first. Its converged value is applied to solve 
the system for the velocity components. Further, since the three linear systems for 
the segregated velocity components are independent, they are solved concurrently 
at each iteration step of CGHS algorithm. After scaling, the number of iterations 
required for convergence of CGHS algorithm is usually a minute fraction of the size 
of the linear system. For example, on a 41^ velocity mesh the CGHS algorithm for 
the linear system for temperature converged in 29 to 41 iterations depending on the 
value of Ra. On the same mesh, the CGHS algorithm for the hnear system for the 
segregated velocity components converged in 37 iterations. 

§9 A numerical Example 

The numerical techniques presented in this report have been applied to simulate nat- 
ural convection in a cubical enclosure containing air with prandtl number 0.71 at Ra 
in the range 10^ to 10^. Initially, the enclosed fluid is stationary and the uniform tem- 
perature in the enclosure and its boundaries is Tc °C. Later, the surface F^ is heated 
uniformly to a temperature °C while the temperature on the surface F^ is held 
fixed at Tc °C. These two surfaces are maintained at these temperatures thereafter 
while the remaining four surfaces of the cube are considered to be perfectly thermally 
insulated. The resulting density variation within the confined fluid is assumed to be 
small enough that the Boussinesq approximation is valid. 
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§10 Results and Discussion 



Steady state results were obtained for the numerical example on unstructured meshes 
with 41^ and 45^ velocity nodes. Steady state solution is assumed when 

< (31) 

where Ui is one of the velocity components and e is taken as 10~^. 
Mesh of up to 41'^ velocity nodes was used for Ra = 10^ and 10^ with At = 1/4000. 
For Ra = 10^ and 10^, when the boundary layer is relatively thinner, non- uniform 
mesh of 45^ velocity nodes was used with At = 1/9000. The time increment in the 
re-discretization of the transport sub-problems was taken in each case as r = 
The construction on the non-uniform mesh, where 13/16 is taken as estimate for the 
boundary layer thickness, is summarized in Table CI 



Ra = 10^ and 10^ 


sub-interval 


No of divisions (coarse mesh) 


45'^ fine mesh 


0<x,y, z< 3/16 


5 




3/16 <x,y, z< 13/16 


12 




13/16 <x,y,z<l 


5 



Table 1: Construction of non-uniform meshes in cubical enclosure. 



Computations were performed on DEC Alpha PW500au, a single processor, virtual 
memory machine and a linux desk-top with 512MB core memory. Computations for 
each Ra were started from the initial conditions given in equationijOj). A mesh with 
41^ velocity nodes requires TOMB of memory while a mesh with 45^ velocity nodes 
requires 99MB of memory. An iteration in time, consisting of solution of the pressure, 
transport and diffusion subproblems, takes average of 1.64 minutes of CPU time. The 
PCG algorithm converged in 7 to 8 iterations after initial transients. 
The Nusselt number, a measure of the dimensionless heat transfer rate across the 
isothermal walls, was computed on the hot wall {x = 0), in terms of the overall 
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Nusselt number, Nu and y-averaged Nusselt number, Nuav{z), defined by the following 
equations: 



For the cubical enclosure Ay = 1.0 and = 1.0 The local heat flux 89/ dx, was 
approximated by a second order forward- difference formula and the integrals were 
evaluated using the trapezoidal rule. 

Variations of all variables with respect to z were investigated at each of the Rayleigh 
numbers applied. These variations in the z-direction are however weaker in magnitude 
than in other directions. In Figure IHl and Figured the distributions of iuy)^^^ {z) at 
X = 0.5 and the y-averaged Nusselt number, Nuav{z) , are illustrated at each Rayleigh 
number. 

§11 Validation 

A comparison of the values of Nu and the mean Nusselt number Nuav{'^-'o) with 
results of Fusegi et al. |^ at each of the Rayleigh numbers used is given in Table El 
The results obtained using this present numerical procedure are in good agreement 
with those of Fusegi et al. Also, in agreement with Fusegi et al. jl], Janssen et al. j3] 
and Mallinson et al. 1], at Ra = 10^ and 10^, the 2;- variations of Nuav{z) are apparent 
near the end walls (2 = and z = 1), where it increases sharply. Janssen et al. [IH] 
also reported a value of 0.2585 for NuRa~^^^ at Ra = 10^. That is, Nu is 8.6396 
at this Rayleigh number and this represents a difference of only 0.27% of the value 
obtained using this present numerical procedure (See Table Ej). 




(32) 




(33) 
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Ka 


Quantity 


Wave Equation 


Fusegi et a/.^ 


% Error 


iU 


I\ u 


i.z400 


i.Uoo 


iz.yb /o 






i.ZOOO 


i.iUo 


10 0/10/ 


iu 


I\ u 


i.y (it 


z.iU 


- O.470 








Z.oUZ 


- (.ZD/o 


10^ 


Nu 


4.2055 


4.361 


- 3.70% 




NUav{0.5) 


4.497 


4.464 


0.73% 


10^ 


Nu 


8.6628 


8.770 


- 1.24% 




NUav{0.5) 


8.8434 


9.012 


- 1.91 % 



Table 2: Comparison of Nusselt numbers with results of Fusegi et al. 



Ra 


Quantity 


Wave Equation 


Janssen et al.^ 


% Error 


10^ 


NuRa^'^/'^ 
Nu 


8.6628 


0.2585 
8.6396 


0.27% 



Table 3: Comparison of Nusselt numbers with results of Janssen et al. 



In agreement with Janssen et al. 0, at Ra = 10^ and 10^, (uy)^^^ (z) has two sharp 
peaks close to the lateral walls. At Ra = 10^, the peaks occur at (0.075, 0.5, 0.15) 
and (0.075, 0.5, 0.85).At Ra = 10^ the peaks occur at (0.0375, 0.5, 0.0937) and 
(0.0375, 0.5, 0.9625). The value and location of the 2;— direction local maximum in 
each case are given in Table IU 





z-direction local maximum 


Location at y = 0.5 


Ra = 10^ 


89.442 


X = 0.075, z = 0.15 


Ra = 10^ 


314.1738 


X = 0.0375, z = 0.0937 



Table 4: Comparison of the peaks for (uy)^^^ (z) with results of Janssen et al. 
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The transverse variations of characteristic quantities, were further investigated. From 
the mesh plots of velocity components, temperature and pressure on the yz-plane at 
X = 0.5 (not shown, see ref. it was observed that for a generic variable /, 

representing u^, Uy, 6 or p, the following relation holds: 

/(x, y, z) ^ /(x, y, 1 - 2;) V X, y, z. (34) 

Uz on the other hand satisfies the relation, 

Uz{x, y, z) ^ -Uz{x, y,l- z) V x, y, z. (35) 

This symmetry property is sometimes imposed by researchers in order to reduce 
calculations to only half of the entire computational domain. 

Also, from the contour plots on xy-plane (at z = 0.5) and mesh plots on x2;-plane 
(at y = 0.5) (ref. 0), it was observed that there is also a special type of symmetry 
about the line (0.5, 0.5, z) (that is, the line y = 0.5 on yz-pla.ne at x = 0.5), through 
the center of gravity, so that the following relations are satisfied: 

u^{x, y, z) ^ -u.j:{1- X, 1-y, z), (36) 

Uy{x,y, z)^ -Uy{l-x,l-y, z), (37) 

Uz{x, y, z) ^ + Uz{l - X, 1-y, z), (38) 

p{x, y, z) ^ +p{l- X, 1-y, z), (39) 

e{x,y,z)^l-e{l-x,l-y,z). (40) 

These two spatial symmetries have been observed by Janssen et al. and they 
exploited it by performing actual computations over only a quarter of the entire 
cubical enclosure. 
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Appendix 



A A Systematic Method of Discretizing A Three-Dimensional Domain 
into Four-Node Tetrahedral Elements 

Let J2h of nodes on Qh and respectively. The first step 

in the domain discretization is the construction of which is illustrated by the 
following steps: 

• Divide the three dimensional domain into brick-like macros, 

1^ {(x, z)\xj < X < Xj + Ax ; Uj < y < yj + Ay ; Zj < z < Zj + Az} , (41) 
where A^2/i = card (^2h)- 

A typical brick-like macro, with corners labeled A, B, C, D, E, F, G, H, is shown 
in Figure El Point A has coordinate {xj, yj, zj). 

• Divide each brick-like macro into two prisms - one with corners H, E, B, C, D, 
A and the other prism with corners H, E, F, B, C, G along diagonals {HC, EB} 
or {DG, AF} (see Figure EI). 

• Sub-divide the first prism into three tetrahedra, through line segments {DE, 
DB, HB}. 

{top vertex-to-bottom edge construction) 

• Sub-divide the second prism into three tetrahedra, through line segments {FH, 
FC, HB}. 

{bottom edge-to-top vertex construction) 

Thus, each 8-cornered brick-like macro on the coarser mesh is divided into six tetra- 
hedral elements for pressure. 

20 



Observe that the three tetrahedral elements on the first prism may be obtained from 
those on second prism by reflecting first prism about the line segment, HC. Thus, we 
may reverse the order of connecting points between these two prisms. That is, we may 
sub-divide the first prism into three tetrahedra through line segments {AH, AC, HB} 
{bottom edge-to-top vertex construction) . The second prism must correspondingly be 
sub-divided into three tetrahedra through line segments {GE, GB, HB} {top vertex- 
to-bottom edge construction) . We may also use diagonal EC in place of HB. If this 
reflection property is taken into consideration, any appropriate combination of line 
segments, as described above, will produce the identical set of six tetrahedral elements. 
This method of discretizing a brick-like macro into six tetrahedra is unique in the sense 
that, the six tetrahedra are generic elements for all the tetrahedral elements in 7^/^ 
and 7/j. The six generic tetrahedra are shown in Figure El 




^ — /'^ 

E F 



Figure 2: A typical brick-like macro. 

The finite set Th for velocity is next constructed from by connecting the edge- 
midpoints on the faces of each tetrahedron and on any resulting vertical mid-plane 
following the steps for constructing the six generic tetrahedra itemized above. 
Thus, each pressure tetrahedral element is sub-divided into eight smaller tetrahedra, 
each of which is similar in shape to one of the six generic tetrahedra depicted in 
Figure IH 

Typical set of 8 velocity elements in two generic pressure macro-elements are delin- 
eated in Figure El 
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gure 3: A systematic division of a brick-like macro into six tetrahedra. 



D DC 




F F 



Figure 4: The six generic tetrahedral elements. 



D ^ C 




Figure 5: Typical pressure macro-elements with 8 sub-tetrahedral. 
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Figure 6: Distribution of (uy) along z-axis in the cubical enclosure at 

(a) Ra = 10^ (b) Ra = 10^ (c) Ra = 10^ and (d) Ra = 10^ 
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Figure 7: Distribution of y-averaged Nusselt number along z-axis in the cubical en- 
closure at (a) Ra = 10^ (b) Ra = 10^ (c) Ra = 10^ and (d) Ra = 10^ 
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